Uma Abordagem Bayesiana com Priore Spike and Slab
Instituto de Matemática, Estatística e Computação Científica (IMECC)
Universidade Estadual de Campinas (UNICAMP)
Considere o modelo: \[ y_i = f(x_i) + e_i \;,\spc i = 1,\dots,n \] com \(e_i \overset{iid}{\sim} \norm(0, \sigma^2)\).
Em Problemas de regressão paramérica, assume-se uma forma funcional para \(f\), por exemplo, \(f(x) = \beta_0 + \beta_1 x + \dots + \beta_p x^p\).
Já em problemas de regressão não paramétrica, \(f\) é desconhecida. Logo, para estimar \(f\) utiliza-se bases como ondaletas, splines entre outras.
Ademais, em termos vetoriais, o modelo é dado por: \[ \bs y = \bs f + \bs e \] onde
\[ \bs y = \begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix}, \hspace{2.5cm} \bs f = \begin{bmatrix}f(x_1)\\\vdots\\f(x_n)\end{bmatrix}, \hspace{2.5cm} \bs e =\begin{bmatrix}e_1\\\vdots\\e_n\end{bmatrix} \]
Vantagens:
Desvantagens:
Seja \(\psi \in \mathbb{L}_2(\RR) = \{f: \RR \to \RR \mid \int f^2 < \infty\}\) uma função que satisfaça a condição \[ C_\psi = \int_\RR \frac{\lvert\Psi(\omega)\rvert^2}{\lvert\omega\rvert} d\omega < \infty \] onde \(\Psi\) é a transformada de Fourier de \(\psi\). Então, dizemos que \(\psi\) é uma ondaleta (mãe).
Além disso, uma vez que a ondaleta foi definida, pode-se gerar ondaletas através de operações de dilatação e translação: \[ \psi_{jk}(x) = 2^{j/2} \psi (2^jx - k) \;, \spc j,k \in \ZZ \]
Como \(\{\psi_{jk}\}_{j,k \in \ZZ}\) formam uma base ortonormal para \(\mathbb{L}_2(\RR)\), qualquer função \(f \in \mathbb{L}_2(\RR)\) pode ser representada por: \[ f(x) = \sum_{j,k \in \ZZ} d_{jk} \psi_{jk}(x) \] onde \(d_{jk}\) é chamado de coeficiente de ondaleta (“detalhes”).
Além disso, também existe a ondaleta pai ou função escala \(\phi\), a qual fazemos translações e dilatações: \[ \phi_{jk}(x) = 2^{j/2} \phi (2^jx - k) \;, \spc j,k \in \NN \]
Com isso, para determinano nível \(j_0 \in \ZZ\) (resolução primária), uma função \(f\) pode ser escrita por: \[ f(x) = \sum_{k \in \ZZ} c_{j_0k} \phi_{j_0k}(x) + \sum_{j \geq j_0} \sum_{k \in \ZZ} d_{jk}\psi_{jk}(x) \] onde \(c_{jk}\) é chamado de coeficiente da função escala.
Do ponto de vista matricial, considerando um vetor de dados \(\bs y \in \RR^n\) diádico, isto é, \(n = 2^J\), \(J \in \NN\). É possível aplicar a transformada discreta de ondaletas (DWT), representada pela matriz \(W\), de forma que \[ \bs d = W \bs y \] onde \(\bs d\) representa os coeficientes empíricos de ondaleta.
Além disso, como \(W\) é uma matriz ortogonal (\(W' = W^{-1}\)), a transformada discreta de ondaletas inversa (IDWT) é dada por: \[ \bs y = W'\bs d \]
Ademais, pela ortogonalidade de \(W\), temos que \(\|\bs d\|^2_2 = \|\bs y\|^2_2\) (Relação de Parseval).
Ondaleta (mãe) de Haar: \[ \psi(x) = \begin{cases} 1 &, \; x \in \left[0, \frac{1}{2}\right)\\ -1 &, \; x \in \left[\frac{1}{2}, 1\right)\\ 0 &, \; c{.}c. \end{cases} \]
Função escala de Haar:
\[ \phi(x) = \begin{cases} 1 &, \; x \in [0, 1]\\ 0 &, \; c{.}c. \end{cases} \]
Agora que fixamos \(\psi\) e \(\phi\), falta apenas definir como calcular \(c_{jk}\) e \(d_{jk}\). Neste caso (considerando a ondaleta de Haar), eles podem ser calculados na primeira aplicação por:
\[ d_{jk} = \frac{1}{\sqrt{2}}(y_{2k} - y_{2k-1}) \hspace{2cm} \text{ e } \hspace{2cm} c_{jk} = \frac{1}{\sqrt{2}}(y_{2k} + y_{2k-1}) \] e, nas demais, por: \[ d_{jk} = \frac{1}{\sqrt{2}}(c_{j+1,2k} - c_{j+1, 2k-1}) \hspace{2cm} \text{ e } \hspace{2cm} c_{jk} = \frac{1}{\sqrt{2}}(c_{j+1,2k} + c_{j+1,2k-1}) \]
Considerando o vetor de dados \(\bs y = (1, 1, 7, 9, 2, 8, 8, 6)'\), a Figura 1 apresenta graficamente a DWT.
Figura 1: Exemplo visual da DWT (Nason, 2008).
A matriz \(W\) que faz essa transformação para \(n=8\) observações é dada por:
\[ W = \begin{bmatrix} \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4}\\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0\\ 0 & 0 &0 & 0 & \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0\\ 0 & 0 & 0 & 0 &0 & 0 & \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}}\\ \frac{-1}{2} & \frac{-1}{2} & \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 &\frac{-1}{2} & \frac{-1}{2} & \frac{1}{2} & \frac{1}{2}\\ \frac{-\sqrt2}{4} & \frac{-\sqrt2}{4} & \frac{-\sqrt2}{4} & \frac{-\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} & \frac{\sqrt2}{4} \end{bmatrix} \]
Portanto, retomando ao modelo de regressão não paramétrica, temos que: \[ \bs y = \bs f + \bs e \] com \(e_i \overset{iid}{\sim} \norm(0, \sigma^2)\).
Então, multiplicando por \(W\) vamos ter o nosso modelo no domínio das ondaletas, dado por: \[ \bs d = \bs\theta + \bs\varepsilon \] onde
Logo, o problema de estimar a função \(f\) foi substituido por estimar os coeficientes de ondaletas \(\bs\theta\) através dos coeficientes empíricos \(\bs d\).
Para isso, a abordagem mais tradicional é a utilização da limiarização (thresholding), a qual consiste em anular coeficientes suficientemente pequenos. As duas propostas mais comuns são a do limiar duro (hard threshold) e limiar suave (soft threshold).
A intuição desta técninca consiste em que o vetor de coficientes de ondaleta normalmente são esparsos. Além disso, coeficientes grandes estão associados às características importantes da função, como mínimos e máximos locais, descontinuidades entre outras. Enquanto que a parte mais suave da função se associa aos coeficientes nulos.
Regra do limiar duro: \[ \delta^H(d) = \begin{cases} 0 &, \text{ se } \lvert d \rvert \leq \lambda\\ d &, \text{ se } \lvert d \rvert > \lambda \end{cases} \]
Figura 1: Regra do limiar duro.
Regra do limiar suave: \[ \delta^S(d) = \begin{cases} 0 &, \text{ se } \lvert d \rvert \leq \lambda\\ \operatorname{sgn}(d) (\lvert d \rvert - \lambda) &, \text{ se } \lvert d \rvert > \lambda \end{cases} \]
Figura 1: Regra do limiar suave.
Existem diversas estratégias para determinar \(\lambda\), algumas delas são:
Validação cruzada.
Universal Threshold: utilizando o MAD (median absolute deviation) para estimar \(\sigma\), temos: \[ \lambda = \sigma \sqrt{2 \ln(n)} \]
SURE Threshold: escolhe o \(\lambda\) que minimiza o estimador não viesado do risco de Stein \(SURE(\lambda, \bs y)\), dado por: \[ SURE(\lambda; \mathbf{y}) = n - 2\#\{i:|y_i| \leq \lambda\} + \sum_{i}(\min\{|y_i|, \lambda\})^2 \]
Figura 1: Resumo do procedimento de estimação.
A razão sinal-ruído (SNR) é dada por:
\[ SNR = \frac{\operatorname{SD}(\text{sinal})}{\operatorname{SD}(\text{ruído})} = \frac{\operatorname{SD}(f)}{\operatorname{SD}(\varepsilon)} \] Para estudos simulacionais, é interessante fixar a razão sinal-ruído para determinar o desepenho da téninca para aquela \(SNR\) e permitir a comparação com outros métodos.
Existem também propostas bayesianas para o problema, por exemplo, a apresentada em Sousa (2020), cuja priori atribuída é: \[\pi(\theta) = \alpha \delta_0(\theta) + (1 - \alpha) g(\theta; \tau)\] onde \(\alpha \in (0,1)\), \(\delta_0\) é o delta de Dirac com massa em \(0\) e \(g(d; \tau)\) é a função densidade de probabilidade logística, dada por: \[ g(\theta; \tau) = \frac{\exp\left\{-\dfrac{\theta}{\tau}\right\}}{\tau \left( 1 + \exp\left\{-\dfrac{\theta}{\tau}\right\} \right)^2} \;\ind_\RR(\theta) \;, \spc\spc \tau > 0 \] Contudo, diferente das estratégias apresentadas anteriormente, propostas bayesianas não zeram os coeficientes, apenas encolhem.
Com isso, é possível estabelecer a distribuição a posteriori \(\theta \mid d\) de cada coeficiente:
\[\begin{align*} \pi(\theta \mid d) &\propto \pi(\theta) \LL(\theta; d)\\ &\propto \left[\alpha \delta_0(\theta) + (1 - \alpha) g(\theta; \tau) \right]\exp\left\{\frac{-1}{2\sigma^2} (d - \theta)^2\right\} \end{align*}\]
A regra de encolhimento é dada pela esperança da posteriori, calculando obtém-se que: \[ \delta(d) = \frac{(1-\alpha)\int_\RR (\sigma u + d) g(\sigma u + d; \tau) \phi(u) du}{\frac{\alpha}{\sigma}\phi\left(\frac{d}{\sigma}\right) + (1 + \alpha) \int_\RR g(\sigma u + d; \tau) \phi(u) du} \] onde \(\phi\) é a densidade da normal padrão, \(u = \dfrac{\theta - d}{\sigma}\) e é utilziado o MAD para estimar \(\sigma\).
Como a integral não é tratável analiticamente, serão utilizados métodos monte carlo para aproximá-la.
A epilepsia é um distúrbio neurológico caracterizado pela ocorrência de convulsões, causadas por uma atividade neuronal excessiva ou sincronizada no cérebro.
Essas crises podem resultar em problemas psiquiátricos, quedas, déficits cognitivos e aumento do risco de morte.
Uma das estratégias mais promissoras é a detecção do estado pré-ictal, período que antecede uma crise.
Identificá-lo pode permitir ações preventivas e reduzir os impactos negativos das convulsões.
Este estudo utiliza dados do eletroencefalograma (EEG) do couro cabeludo de \(14\) pacientes, registrados com eletrodos reutilizáveis de prata e ouro.
Os dados foram retirados de Detti et al. (2020), da Unidade de Neurologia e Neurofisiologia da Universidade de Siena, Itália.
O banco apresenta variáveis como gênero, idade, número de canais de EEG, quantidade de crises e o tempo das gravações (em minutos).
Para este estudo, foi seleiconado o paciente \(5\) com tamanho amostral correspondente a \(n = 2^{15} = 32.768\) observações.
Nason, G.P. (2008). Wavelet Methods in Statistics with R. Springer.
Nason, G.P. (2024). wavethresh: Wavelets Statistics and Transforms. R package version 4.7.3, https://CRAN.R-project.org/package=wavethresh.
Sousa, A.R.S. (2018). Encolhimento Bayesiano Sob Priori Beta de Coeficientes de Ondaletas em Modelos com Erros Gaussianos e Positivos. Tese (doutorado) - Universidade Estadual de Campinas, Instituto de Matemática Estatística e Computação Científica, Campinas, SP. DOI: 10.47749/T/UNICAMP.2018.1080773.
Sousa, A.R.S. (2020). Bayesian wavelet shrinkage with logistic prior. Communications in Statistics - Simulation and Computation. DOI: 10.1080/03610918.2020.1747076.
Detti, P., Vatti, G., Lara, Z.M. (2020). EEG Synchronization Analysis for Seizure Prediction: A Study on Data of Noninvasive Recordings. Processes, 8(7), 846. DOI: 10.3390/pr8070846.